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Jeffrey Tsui and Sam Reisenfeld 


Abstract— A highly accurate DFT-based complex exponential 
parameter estimation algorithm is presented in this paper. It 
will be shown that for large number of samples and high sig- 
nal to noise ratio (SNR), the phase estimation error variance 
performance is only 0.0475 dB above the Cramer-Rao lower 
bound (CRLB) for phase estimation with unknown frequency 
and phase. The amplitude estimation error variance perfor- 
mance was found to lay on the CRLB for amplitude estima- 
tion. Exact phase and amplitude estimation can be achieved 
in the noiseless case with this algorithm. The algorithm has 
low implementation computational complexity and is suitable 
for numerous real time digital signal processing applications. 


Keywords— frequency estimation, phase estimation, amplitude 
estimation, DFT-based parameter estimation, spectral estima- 
tion, digital signal processing algorithm, complex exponential 
parameter estimation. 


1. Introduction 


Frequency, phase and amplitude estimation of a complex 
exponential is classical problem in statistical signal pro- 
cessing. The precision of the phase and amplitude estimate 
is directly related to the accuracy of the frequency estimate. 
There are two major classes of complex exponential fre- 
quency estimation algorithm in the existing literature. The 
first class is the classical discrete Fourier transform (DFT) 
and phase averager based frequency estimation algorithms 
such as [3] and [4]. This class of estimation algorithm 
is very computationally efficient, however they suffer from 
poor error performance at low signal to noise ratio (SNR). 
The second class is the parametric frequency estimation 
algorithms such as the very popular MUSIC and ESPRIT 
algorithms [5, 6]. This class of estimation algorithm has 
very good error performance across a wider range of SNR, 
however, they are very computationally intensive and not 
suitable for real-time application. Obviously, it would be 
ideal if there exists a complex exponential parameter esti- 
mator that is both computationally efficient and has good 
performance at low SNR. 

In 2003, Reisenfeld and Aboutanios [2] discovered a fre- 
quency estimator that can precisely estimating the fre- 
quency of a complex exponential in an iterative manner 
by applying contraction mapping method on two modified 
DFT coefficients. In a subsequence publication in 2004, 
Reisenfeld [1] further enhanced the previous published al- 
gorithm by improving the convergence property of the es- 
timator. It can be shown that this particular frequency es- 
timator can yield a frequency estimate that is 0.0633 dB 
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of the Cramer-Rao lover band (CRLB) with approximately 
Nlog,N +4N complex multiplications, where N is the 
number of samples used representing the signal. 

In this paper, it will be shown that combining the said 
frequency estimator with maximum likelihood (ML) phase 
and amplitude estimators yields a highly accurate param- 
eter estimator for complex exponentials. In the noiseless 
case, it is possible to obtain the exact phase and amplitude 
estimates with this estimator. In additive white Gaussian 
noise (AWGN) channel, the said estimates approach very 
close the their respective CRLB at relative high SNR. The 
relationship between the number of samples, N, and the 
operating point of the parameter estimator in terms of SNR 
will be given. It will also be shown that it is possible to 
use the said relationship to further optimise the computa- 
tion complexity of the estimator. 

The rest of the paper will be organised as follows: Sec- 
tion 2 introduces the proposed DFT-based parameter esti- 
mator, Section 3 will provide the performance analysis of 
the proposed parameter estimator. Section 4 will describe 
further enhancements that can be made to the original fre- 
quency estimator in [1] to reduce its computation complex- 
ity. Section 5 contains the simulated results of the perfor- 
mance of the proposed algorithm and this will be followed 
by the conclusion. 


2. The DFT-based parameter estimator 


The proposed DFT-based parameter estimator involves 
a two stage process. First, the frequency of the received 
carrier is estimated by the frequency estimator as described 
in [1]. The frequency estimate obtained is then used to 
eliminate the frequency component of the carrier, leaving 
only the phase and amplitude component to be estimated 
by the ML phase and amplitude estimator. 


2.1. The DFT-based complex exponential frequency 
estimator 


Consider a complex exponential, r[n], with amplitude, A;, 
frequency, f- € [0, f;), and phase 0. € [0,27). Mathemati- 
cally, r[n] can be represented as 


r[n] = Ace #4 %) 4. 1 In], (1) 


where n= 0,1,2,.....N—1, T; = 1/fs is the sampling pe- 
riod, and 7 [n] is a sequence of independent complex Gaus- 
sian variables with mean of zero and variance o7. 
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In noiseless condition, the magnitude spectra the DFT of 
Eq. (1) will have even symmetry about its frequency. The 
recursive algorithm as described in [1] exploits this condi- 
tion by employing a discriminate that works on a contrac- 
tion principle in minimizing the difference in the magnitude 
of two modified DFT coefficients which are plus and minus 
half a DFT bin away from an estimated frequency. Due to 
the even symmetric nature of the magnitude spectra, the 
difference in the magnitude of the two modified DFT coef- 
ficients will eventually reduced to zero in the noiseless case 
as the estimated frequency approaches the true frequency 
with the increasing number of recursions of the algorithm. 
The frequency can then be estimated as the mean of the 
frequencies of the modified DFT coefficients at which dif- 
ference in their magnitude equals to zero. 

Summarizing the DFT-based frequency estimator in [1]: 


1. Perform a coarse frequency estimate such as the 
one described in Rife and Boorstyn [3], in which 
{r[n] y is the input to a N point complex DFT 
and a peak search is done on the magnitudes of 
the DFT output coefficients, to obtain the initial fre- 
quency estimate, fo. This estimate is obtained by, 
fo =kmaxfs ‘ N, where kyax is the index of the maxi- 
mum magnitude DFT output coefficient. 


2. Calculate the modified DFT coefficients Q%, and Bn 
defined by 


N-1 oo. F 
Om = Yo r[nje Perna), (2) 
n=0 


N-1 a ee 
Bn = y r|n] e Pan (Sn Ts+ ay), (3) 


n=0 


3. Calculate the discriminate D,, defined as 


in| — | Om 
Dm = |Bm| | nl (4) 


7 [Brn + |Qn| 


4. Calculate the new adjusted frequency with the for- 
mula: 


1 


Frat = fnt+ tan [Pmtan (=) | fs. (5) 


5. Perform Steps 2-4 recursively for m= 1, 2, 3, ..., 
M-1. 


2.2. Combining the frequency estimator with ML phase 
and amplitude estimator 


Using the frequency estimate f, obtained from the fre- 
quency estimator described above, it is possible to elimi- 
nate the frequency component of r [n] by multiplying it with 
the conjugate of the complex exponential with a frequency 
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of fm. Denoting r4.g[n] as the result of the multiplication 
we have 


TA. [n| = (Ach@nsenti+ée) + n inl) @i2tfmnTs 
= Ael(2enTs +6) + n [n] ; (6) 


where € = f. = 

It was shown in [7] that the ML phase estimate of a complex 

exponential is obtained by taking its arctangent. Taking the 
N 


-1 
arctangent of Y r4.9 [n] yields 
n=0 


n=0 


im [E (astomre +n) 
-1 


Re 's (Aei(2nenT, +6.) +n [n]) 


N-1 
¥ Asin (2menT, + 0.) + 1s 
n=0 


py «MD 


N-1 
¥ Acos (27enT; + 8) + Ne 
n=0 


where 


N-1 N 
Nc =Re (x 7 ni) , E[n-] =0, Var|ne] = 9) 
n=0 


N-1 N 2 
n=im(¥ ntl) E[n]=0, VarinJ=——, ©) 
n=0 


and 6.. fm 1S the phase estimate based upon the estimated 
frequency. 


It was also shown in [7] that the ML amplitude estimate 
of a complex exponential is obtained by taking its absolute 


N-1 
value. Taking the absolute value of Yr, 9 [n] yields 
n=0 


2 
+Im 


2 


Yi (Aecienents+6) +0 [n)) 
n=0 


1 
Aj, = 7 = [s(acernan (nl) 


1=0 


2 2 


N-1 
+] } Asin (2menT,) +n, 


n=0 


1 | Pye 
= i p Acos (27enT,) + Nc 


n=0 


(10) 


Notice the amplitude estimation is not affected by the phase 
angle of the complex exponential 0,. 

From Eqs. (7) and (10), it is obvious that a highly accurate 
frequency estimator can assist in reducing the error and 
improve on the accuracy of the phase estimation. Hence 
the described frequency estimator is very suitable for this 
joint estimation of phase and amplitude due to its superior 
error variance performance. 
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3. Performance of the parameter 
estimator 
3.1. Frequency estimator 
The author in [1] has proven this algorithm only requires 
two iterations for the variance of the frequency estimate / 


to converge to less than or equal to 0.063 dB above the 
CRLB for frequency estimation. 


3.2. Phase estimator 


Assuming high SNR, using Taylor series expansion as de- 
scribed in [8] the variance of Eq. (7) was found to be 


i N2(N—1) sin? (4) tan? (4) +2 
var |6, ;, _ ( ) (sy) C3) ; (11) 
C,Jm 4oN 
where p, the SNR equals to 
Az 
p= a (12) 


Since this algorithm jointly estimates the frequency and 
phase of the observed signal, it is appropriate to compare 
the performance of the phase estimator to the CRLB of 
phase estimation with unknown frequency and phase which 
is given by [7] 


(2N —1) 


~ Np(Waly “ 


CRLB joint (8) 


Therefore 


Var a 
CRLB joint (0) 


(1? (W— 1)? sin? (4) tan? (5) +2) (N+1) 
> 4(2N — 1) =) 


For large p and large N, it can be shown, 


Var [4.4] 
lim 101 CRIB on (6) | 
Piss 0810 CRLB joint (0) 


4 
1 
10log1o (4+3) = 0.0475 dB. (15) 


Figure | shows the convergence property of Var [ @ese, f| to 
the CRLB as a function of the number of samples N. It can 
be seen that variance of the phase estimator deviates from 
the CRLB as the number of samples N increases and ap- 
proaches the asymptotic limit of 0.0475 dB. This due to the 
convergence behavior of the frequency estimator as stated 
in [1] where for small N, less information is discarded for 
not using all the DFT coefficients hence the performance 
degradation is less than when N is large. 
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Var |, jn] 


Fig. 1. CREB (8) 


in dB as a function of number of samples, NV. 


3.3. Amplitude estimator 


Assuming high SNR, using Taylor series expansion as de- 
scribed in [8] the variance of Eq. (10) was found to be 


; o 
var |A;,| = (16) 


which agrees with the CRLB derived by Rife and Boorstyn 
in [3]. 


4. Further enhancements 


As describe in Subsection 2.1, the frequency estimator 
in [1] requires an initial frequency estimate obtained by 
performing a fast Fourier transform (FFT) operation on the 
signal samples. Since FFT requires Nlog, N complex mul- 
tiplications, it is desirable to keep the number of samples NV 
as low as possible for the initial coarse estimate. How- 
ever, the CRLB for frequency estimation monotonically de- 
creases with the increasing number of samples used for the 
estimation at a fixed SNR. Hence, reducing the number of 
samples used for the frequency estimate may not produce 
an estimate that will meet the required accuracy. 

One way to reduce the complexity of the estimator and 
obtain the required accuracy is to modify the frequency 
estimator so that it beings the initial coarse frequency 
estimation with a low number of samples and dynamically 
increase the number of samples for each iteration of fre- 
quency estimation algorithm. The modified version of the 
DFT-based frequency estimator is summarised as follows: 


Denoting N;, where i = 0,1,2,3,...,7, as the number of 
samples used in the ith pass of the frequency estimation 
algorithm. Note that the number of samples, N;, for each 
pass must satisfy the following relationship No < Ni < N2 
<... <M, and N; = N which is the number of samples 
required to obtain the desire accuracy. 


1. Perform a coarse frequency estimate such as the 
one described in Rife and Boorstyn [3], in which 
{r[n] at is the input to a No point complex DFT 
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and a peak search is done on the magnitudes of 
the DFT output coefficients, to obtain the initial fre- 
quency estimate, fo. This estimate is obtained by, 
fo =kmaxfs jf No, Where kmax 1s the index of the max- 
imum magnitude DFT output coefficient. 


2. Perform the ith pass of the frequency estimator which 
consists of the following steps: 
2.1. Recursion is started at m= 0. 
2.2. Set N; = No. 


2.3. Calculate the modified DFT coefficients Qj, 
and f,, defined by 


NG oe) -i2an( fnts— ah) 

Om = ) r{nje se a a (17) 
n=0 
= —j2an (fn T+ ! ) 

Bn = ¥- r(nje ea) (18) 
n=0 


2.4. Calculate the discriminate D,, defined as 


|Bre| — | On| 
Dyn = = ——. (19) 
" [Bin| + [On| 


2.5. Calculate the new adjusted frequency with the 


formula: 
4 
wn) |f 


(20) 


7 1 _ 
Fin+1 = Tin + = tan 1 >, tan ( 


2.6. Perform Steps 2.3-2.5 recursively for m = 
1,2,3,..M =I, 


3. Set the value of fo to the value of = as found in 
Step 2.5. Repeat Step 2 for Nj = Ni, Ni =No,..., 
N; = Ny. 


Choosing the value of N;’s. Thus far, the discussion has 
been on the reducing the number of samples to save com- 
putational complexity. However, the question of how one 
practically chose the values of N;’s remains. The choices 
of the value of the N;’s are govern by the frequency er- 
ror variance of the estimate from the individual ith passes. 
If the frequency estimate, f,,, of the ith was too far away 
from actual frequency, it will cause the frequency estimator 
converge on an incorrect frequency. In fact, in the original 
algorithm presented in [1], the initial frequency estimate 
error must be bound between the frequency that is repre- 
sented by +1 / 2 of a DFT bin to ensure convergence of the 
algorithm. Since we are increasing the number of samples, 
Ni, at each pass, one must ensure the frequency estimate 
error of the ith pass must be smaller than +1 j 2 the fre- 
quency represented by a DFT bin of the next pass. This 
relationship can be mathematically represented as 


a) 
N;sin (+4; tan? (#) 1 
4p 7? s 2Niv1 en) 
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In addition to the above condition, the value of No, which 
corresponds to the length of initial FFT for the initial peak 
search, has an extra constraint in the form of the perfor- 
mance threshold at low SNR as discussed in [3]. The 
performance threshold has to do with the nonlinear nature 
of the frequency estimation problem as low SNR. A detail 
discussion on the relationship between performance thresh- 
old, SNR and the number of samples is out of the scope 
of this paper. There are a number of papers that addresses 
this issue and readers are recommended to look at refer- 
ences [9-12] for detail analysis of performance threshold. 
Only the findings from [9] are discussed in this paper be- 
cause of the simplistic nature of the results and the ease of 
applying them to determine the optimum value of No given 
it has to operate above certain SNR. 

In [9], the author derived an approximate threshold indica- 
tor given by 


E(N&*)? = ——, (22) 


where 6* is the approximation of the normalised frequency 
error. It was found there is a relationship between the indi- 
cator quantity given in Eq. (22) and the mean square (MS) 
phase error given by 

3 


3x2 

E(NS&*)? = —— =~E(6)°. 23 
(NS")" = 5 om a (6) (23) 

The author in [9] also found the MS phase error associated 


with the threshold is roughly E (6)” = 0.0625 rad2. Table 1 
was built using Eq. (23) at common values of N. 


Table 1 
SNR threshold for various values of N 
N Threshold SNR [dB] 
1024 —15.05 
512 —12.04 
256 —9.03 
128 —6.02 
64 —3.01 
32 0 


Since these threshold values are approximations, one should 
allow a operating margin of at least 1.5 dB when deciding 
upon the value of No. For example, if the requirement is 
for the estimator to operate at —3 dB, one would choose 
No = 128 over No = 64 to guarantee proper operation 
at —3 dB. Compare the values in Table | to the relation- 
ship as stated in Eq. (21), one can conclude the performance 
threshold will dominate in deciding on the value of No. 


5. Simulation results 


Figures 2, 3 and 4 shows the simulated results of the error 
variance performance as a function of SNR for the fre- 
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Fig. 2. The error variance of the DFT-based frequency estimator 
as a function of SNR. 
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Fig. 3. The error variance of the phase estimator as a function 
of SNR. 
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Fig. 4. The error variance of the amplitude estimator as a func- 
tion of SNR. 
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quency, phase and amplitude estimator described in Sec- 
tion 2 compared to results obtained by the Matlab™ “root- 
music” algorithm and TLS-ESPRIT algorithm presented 
in [13]. The number of data points used in the simula- 
tion equals to 256 and the size of the autocorrelation ma- 
trix used for the “rootmusic” and TLS-ESPRIT algorithm 
equals to 64. For each trial a random frequency and phase 
is generated from two independent uniform distributions 
with the range —f; / 2 to fs / 2 and 0 to 27, respectively. 
The results shown are obtained by averaging over 6000 tri- 
als. It can be seen that the simulated results obtained by 
the proposed algorithm is on par with those obtained by us- 
ing “rootmusic” and TLS-ESPRIT and yet computationally 
very efficient. 

Figure 5 shows a comparison of the frequency estimation 
error variance as a function of SNR between the modified 
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Fig. 5. The error variance of the modified frequency estimator 
compared to the original frequency estimator. 


frequency estimator as described in Section 4 against the 
frequency estimator as described in Subsection 2.1. The 
frequency estimate obtained from two passes of the modi- 
fied frequency estimator with No = 128 and N; = 512. The 
number of samples was kept constant at N = 512 for the 
frequency estimator as described in Subsection 2.1. The 
results shown are obtained by averaging over 6000 trials. 
As expected, the modified version of the frequency esti- 
mator’s performance threshold is at a higher SNR than the 
original estimator. However the modified frequency estima- 
tor only required No log, No +4No +4N 1 = 3456 complex 
multiplications which is approximately half of what is re- 
quired for the original frequency estimator which equals to 
Nlog, N+4N = 6656 complex multiplications. 


6. Conclusion 


Expanding upon [1], a new algorithm for joint frequency, 
phase and amplitude estimation of a complex exponential 
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has been presented. It was shown that given sufficient 
number of samples, at high SNR the variance of the phase 
estimator approaches the asymptotic limit of 0.0475 dB 
above the CRLB. It was also shown that the modified ver- 
sion of the frequency estimator can significantly reduce 
the computation complexity with compromising the over- 
all error variance performance except increasing the oper- 
ation SNR threshold. With the advantage of being compu- 
tationally efficient, this type of joint frequency and phase 
estimator is well suited for real time application such as 
timing and carrier synchronization. 


Appendix 


Variance of the DFT-based phase estimator 


The arctangent of ra 9 [n] as defined in Eq. (6) was given 
as 


N-1 
)y) Asin (22enT; + 0c) + 1s 


8.5, jag! — . (24) 
> Acos (2€nT; + 8) + Ne 
n=0 


From [1], € the frequency estimation error is a random 
variable with zero mean and variance of 


Nsin2 (=) tan? (=) 
2N 2N 
422SNR : 


Var [€] = (25) 


Again, the variance of 6. fe in Eq. (24) can be found by 
using the technique of linearization of function of random 
variables presented in [8]. From the structure of the dis- 
criminate D,, defined in Eq. (4) one can conclude € and 1, 
are uncorrelated and € and 1s; are uncorrelated. Expand- 
ing Eq. (24) in a three dimensional Taylor series expansion 
with respect to these variable gives 


~ 2 x 2 
6 > « 
Var | in| = 6 Var le] + on. Var [Ne] 
= 2 
+ | eh | varing] (26) 
Ons dl 
where 
06 >; 06 ; 
Cif 2 cfm _ _ 
de | \ Oe et alae 
e=E[e] ,Ne=E[Ne] Ns=E[Ns] 
(27) 
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d 6.., _ a 6,2. _ sin (0) 
re) Cc — 0 is aa TA ’ 
" q e=Ele|, Ne=E [Ne], Ns=E [Ns] - 
(28) 
J i = 8.3, = cos (8) 
ans Ons NA 
: ™ lee le]. Ne=E[Ne].Ns=E [Ns] 
(29) 
Solving Eq. (26) yields 
. N? (N — 1)? sin? (4) tan? (4) +2 
Var 844] = 4N SNR Soe 
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